home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / numerical / slatec / zbknu.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  24.8 KB  |  692 lines
  1. ;;; Compiled by f2cl version 2.0 beta 2002-05-06
  2. ;;; 
  3. ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
  4. ;;;           (:coerce-assigns :as-needed) (:array-type ':simple-array)
  5. ;;;           (:array-slicing nil) (:declare-common nil)
  6. ;;;           (:float-format double-float))
  7.  
  8. (in-package "SLATEC")
  9.  
  10.  
  11. (let ((kmax 30)
  12.       (czeror 0.0)
  13.       (czeroi 0.0)
  14.       (coner 1.0)
  15.       (conei 0.0)
  16.       (ctwor 2.0)
  17.       (r1 2.0)
  18.       (dpi 3.141592653589793)
  19.       (rthpi 1.2533141373155003)
  20.       (spi 1.909859317102744)
  21.       (hpi 1.5707963267948966)
  22.       (fpi 1.8976999933151775)
  23.       (tth 0.6666666666666667)
  24.       (cc (make-array 8 :element-type 'double-float)))
  25.   (declare (type (simple-array double-float (8)) cc)
  26.            (type double-float tth fpi hpi spi rthpi dpi r1 ctwor conei coner
  27.             czeroi czeror)
  28.            (type f2cl-lib:integer4 kmax))
  29.   (f2cl-lib:fset (f2cl-lib:fref cc (1) ((1 8))) 0.5772156649015329)
  30.   (f2cl-lib:fset (f2cl-lib:fref cc (2) ((1 8))) -0.04200263503409524)
  31.   (f2cl-lib:fset (f2cl-lib:fref cc (3) ((1 8))) -0.04219773455554433)
  32.   (f2cl-lib:fset (f2cl-lib:fref cc (4) ((1 8))) 0.007218943246663099)
  33.   (f2cl-lib:fset (f2cl-lib:fref cc (5) ((1 8))) -2.1524167411495104e-4)
  34.   (f2cl-lib:fset (f2cl-lib:fref cc (6) ((1 8))) -2.013485478078824e-5)
  35.   (f2cl-lib:fset (f2cl-lib:fref cc (7) ((1 8))) 1.1330272319816956e-6)
  36.   (f2cl-lib:fset (f2cl-lib:fref cc (8) ((1 8))) 6.116095104481417e-9)
  37.   (defun zbknu (zr zi fnu kode n yr yi nz tol elim alim)
  38.     (declare (type (simple-array double-float (*)) yr yi)
  39.              (type f2cl-lib:integer4 kode n nz)
  40.              (type double-float zr zi fnu tol elim alim))
  41.     (prog ((cssr (make-array 3 :element-type 'double-float))
  42.            (csrr (make-array 3 :element-type 'double-float))
  43.            (bry (make-array 3 :element-type 'double-float))
  44.            (cyr (make-array 2 :element-type 'double-float))
  45.            (cyi (make-array 2 :element-type 'double-float)) (i 0) (iflag 0)
  46.            (inu 0) (k 0) (kflag 0) (kk 0) (koded 0) (idum 0) (j 0) (ic 0)
  47.            (inub 0) (nw 0) (aa 0.0) (ak 0.0) (ascle 0.0) (a1 0.0) (a2 0.0)
  48.            (bb 0.0) (bk 0.0) (caz 0.0) (cbi 0.0) (cbr 0.0) (cchi 0.0)
  49.            (cchr 0.0) (cki 0.0) (ckr 0.0) (coefi 0.0) (coefr 0.0) (crscr 0.0)
  50.            (csclr 0.0) (cshi 0.0) (cshr 0.0) (csi 0.0) (csr 0.0) (czi 0.0)
  51.            (czr 0.0) (dnu 0.0) (dnu2 0.0) (etest 0.0) (fc 0.0) (fhs 0.0)
  52.            (fi 0.0) (fk 0.0) (fks 0.0) (fmui 0.0) (fmur 0.0) (fr 0.0) (g1 0.0)
  53.            (g2 0.0) (pi_ 0.0) (pr 0.0) (pti 0.0) (ptr 0.0) (p1i 0.0) (p1r 0.0)
  54.            (p2i 0.0) (p2m 0.0) (p2r 0.0) (qi 0.0) (qr 0.0) (rak 0.0) (rcaz 0.0)
  55.            (rzi 0.0) (rzr 0.0) (s 0.0) (smui 0.0) (smur 0.0) (sti 0.0)
  56.            (str 0.0) (s1i 0.0) (s1r 0.0) (s2i 0.0) (s2r 0.0) (tm 0.0) (t1 0.0)
  57.            (t2 0.0) (elm 0.0) (celmr 0.0) (zdr 0.0) (zdi 0.0) (as 0.0)
  58.            (alas 0.0) (helim 0.0))
  59.       (declare (type (simple-array double-float (2)) cyi cyr)
  60.                (type (simple-array double-float (3)) cssr csrr bry)
  61.                (type double-float helim alas as zdi zdr celmr elm t2 t1 tm s2r
  62.                 s2i s1r s1i str sti smur smui s rzr rzi rcaz rak qr qi p2r p2m
  63.                 p2i p1r p1i ptr pti pr pi_ g2 g1 fr fmur fmui fks fk fi fhs fc
  64.                 etest dnu2 dnu czr czi csr csi cshr cshi csclr crscr coefr
  65.                 coefi ckr cki cchr cchi cbr cbi caz bk bb a2 a1 ascle ak aa)
  66.                (type f2cl-lib:integer4 nw inub ic j idum koded kk kflag k inu
  67.                 iflag i))
  68.       (setf caz (zabs zr zi))
  69.       (setf csclr (/ 1.0 tol))
  70.       (setf crscr tol)
  71.       (f2cl-lib:fset (f2cl-lib:fref cssr (1) ((1 3))) csclr)
  72.       (f2cl-lib:fset (f2cl-lib:fref cssr (2) ((1 3))) 1.0)
  73.       (f2cl-lib:fset (f2cl-lib:fref cssr (3) ((1 3))) crscr)
  74.       (f2cl-lib:fset (f2cl-lib:fref csrr (1) ((1 3))) crscr)
  75.       (f2cl-lib:fset (f2cl-lib:fref csrr (2) ((1 3))) 1.0)
  76.       (f2cl-lib:fset (f2cl-lib:fref csrr (3) ((1 3))) csclr)
  77.       (f2cl-lib:fset (f2cl-lib:fref bry (1) ((1 3)))
  78.                      (/ (* 1000.0 (f2cl-lib:d1mach 1)) tol))
  79.       (f2cl-lib:fset (f2cl-lib:fref bry (2) ((1 3)))
  80.                      (/ 1.0 (f2cl-lib:fref bry (1) ((1 3)))))
  81.       (f2cl-lib:fset (f2cl-lib:fref bry (3) ((1 3))) (f2cl-lib:d1mach 2))
  82.       (setf nz 0)
  83.       (setf iflag 0)
  84.       (setf koded kode)
  85.       (setf rcaz (/ 1.0 caz))
  86.       (setf str (* zr rcaz))
  87.       (setf sti (* (- zi) rcaz))
  88.       (setf rzr (* (+ str str) rcaz))
  89.       (setf rzi (* (+ sti sti) rcaz))
  90.       (setf inu (f2cl-lib:int (+ fnu 0.5)))
  91.       (setf dnu (- fnu inu))
  92.       (if (= (abs dnu) 0.5) (go label110))
  93.       (setf dnu2 0.0)
  94.       (if (> (abs dnu) tol) (setf dnu2 (* dnu dnu)))
  95.       (if (> caz r1) (go label110))
  96.       (setf fc 1.0)
  97.       (multiple-value-bind
  98.           (var-0 var-1 var-2 var-3 var-4)
  99.           (zlog rzr rzi smur smui idum)
  100.         (declare (ignore var-0 var-1))
  101.         (setf smur var-2)
  102.         (setf smui var-3)
  103.         (setf idum var-4))
  104.       (setf fmur (* smur dnu))
  105.       (setf fmui (* smui dnu))
  106.       (multiple-value-bind
  107.           (var-0 var-1 var-2 var-3 var-4 var-5)
  108.           (zshch fmur fmui cshr cshi cchr cchi)
  109.         (declare (ignore var-0 var-1))
  110.         (setf cshr var-2)
  111.         (setf cshi var-3)
  112.         (setf cchr var-4)
  113.         (setf cchi var-5))
  114.       (if (= dnu 0.0) (go label10))
  115.       (setf fc (* dnu dpi))
  116.       (setf fc (/ fc (sin fc)))
  117.       (setf smur (/ cshr dnu))
  118.       (setf smui (/ cshi dnu))
  119.      label10
  120.       (setf a2 (+ 1.0 dnu))
  121.       (setf t2
  122.               (exp
  123.                (-
  124.                 (multiple-value-bind
  125.                     (ret-val var-0 var-1)
  126.                     (dgamln a2 idum)
  127.                   (declare (ignore var-0))
  128.                   (setf idum var-1)
  129.                   ret-val))))
  130.       (setf t1 (/ 1.0 (* t2 fc)))
  131.       (if (> (abs dnu) 0.1) (go label40))
  132.       (setf ak 1.0)
  133.       (setf s (f2cl-lib:fref cc (1) ((1 8))))
  134.       (f2cl-lib:fdo (k 2 (f2cl-lib:int-add k 1))
  135.                     ((> k 8) nil)
  136.         (tagbody
  137.           (setf ak (* ak dnu2))
  138.           (setf tm (* (f2cl-lib:fref cc (k) ((1 8))) ak))
  139.           (setf s (+ s tm))
  140.           (if (< (abs tm) tol) (go label30))
  141.          label20))
  142.      label30
  143.       (setf g1 (- s))
  144.       (go label50)
  145.      label40
  146.       (setf g1 (/ (- t1 t2) (+ dnu dnu)))
  147.      label50
  148.       (setf g2 (* (+ t1 t2) 0.5))
  149.       (setf fr (* fc (+ (* cchr g1) (* smur g2))))
  150.       (setf fi (* fc (+ (* cchi g1) (* smui g2))))
  151.       (multiple-value-bind
  152.           (var-0 var-1 var-2 var-3)
  153.           (zexp fmur fmui str sti)
  154.         (declare (ignore var-0 var-1))
  155.         (setf str var-2)
  156.         (setf sti var-3))
  157.       (setf pr (/ (* 0.5 str) t2))
  158.       (setf pi_ (/ (* 0.5 sti) t2))
  159.       (multiple-value-bind
  160.           (var-0 var-1 var-2 var-3 var-4 var-5)
  161.           (zdiv 0.5 0.0 str sti ptr pti)
  162.         (declare (ignore var-0 var-1 var-2 var-3))
  163.         (setf ptr var-4)
  164.         (setf pti var-5))
  165.       (setf qr (/ ptr t1))
  166.       (setf qi (/ pti t1))
  167.       (setf s1r fr)
  168.       (setf s1i fi)
  169.       (setf s2r pr)
  170.       (setf s2i pi_)
  171.       (setf ak 1.0)
  172.       (setf a1 1.0)
  173.       (setf ckr coner)
  174.       (setf cki conei)
  175.       (setf bk (- 1.0 dnu2))
  176.       (if (or (> inu 0) (> n 1)) (go label80))
  177.       (if (< caz tol) (go label70))
  178.       (multiple-value-bind
  179.           (var-0 var-1 var-2 var-3 var-4 var-5)
  180.           (zmlt zr zi zr zi czr czi)
  181.         (declare (ignore var-0 var-1 var-2 var-3))
  182.         (setf czr var-4)
  183.         (setf czi var-5))
  184.       (setf czr (* 0.25 czr))
  185.       (setf czi (* 0.25 czi))
  186.       (setf t1 (* 0.25 caz caz))
  187.      label60
  188.       (setf fr (/ (+ (* fr ak) pr qr) bk))
  189.       (setf fi (/ (+ (* fi ak) pi_ qi) bk))
  190.       (setf str (/ 1.0 (- ak dnu)))
  191.       (setf pr (* pr str))
  192.       (setf pi_ (* pi_ str))
  193.       (setf str (/ 1.0 (+ ak dnu)))
  194.       (setf qr (* qr str))
  195.       (setf qi (* qi str))
  196.       (setf str (- (* ckr czr) (* cki czi)))
  197.       (setf rak (/ 1.0 ak))
  198.       (setf cki (* (+ (* ckr czi) (* cki czr)) rak))
  199.       (setf ckr (* str rak))
  200.       (setf s1r (+ (- (* ckr fr) (* cki fi)) s1r))
  201.       (setf s1i (+ (* ckr fi) (* cki fr) s1i))
  202.       (setf a1 (* a1 t1 rak))
  203.       (setf bk (+ bk ak ak 1.0))
  204.       (setf ak (+ ak 1.0))
  205.       (if (> a1 tol) (go label60))
  206.      label70
  207.       (f2cl-lib:fset (f2cl-lib:fref yr (1) ((1 n))) s1r)
  208.       (f2cl-lib:fset (f2cl-lib:fref yi (1) ((1 n))) s1i)
  209.       (if (= koded 1) (go end_label))
  210.       (multiple-value-bind
  211.           (var-0 var-1 var-2 var-3)
  212.           (zexp zr zi str sti)
  213.         (declare (ignore var-0 var-1))
  214.         (setf str var-2)
  215.         (setf sti var-3))
  216.       (multiple-value-bind
  217.           (var-0 var-1 var-2 var-3 var-4 var-5)
  218.           (zmlt s1r s1i str sti (f2cl-lib:fref yr (1) ((1 n)))
  219.            (f2cl-lib:fref yi (1) ((1 n))))
  220.         (declare (ignore var-0 var-1 var-2 var-3))
  221.         (f2cl-lib:fset (f2cl-lib:fref yr (1) ((1 n))) var-4)
  222.         (f2cl-lib:fset (f2cl-lib:fref yi (1) ((1 n))) var-5))
  223.       (go end_label)
  224.      label80
  225.       (if (< caz tol) (go label100))
  226.       (multiple-value-bind
  227.           (var-0 var-1 var-2 var-3 var-4 var-5)
  228.           (zmlt zr zi zr zi czr czi)
  229.         (declare (ignore var-0 var-1 var-2 var-3))
  230.         (setf czr var-4)
  231.         (setf czi var-5))
  232.       (setf czr (* 0.25 czr))
  233.       (setf czi (* 0.25 czi))
  234.       (setf t1 (* 0.25 caz caz))
  235.      label90
  236.       (setf fr (/ (+ (* fr ak) pr qr) bk))
  237.       (setf fi (/ (+ (* fi ak) pi_ qi) bk))
  238.       (setf str (/ 1.0 (- ak dnu)))
  239.       (setf pr (* pr str))
  240.       (setf pi_ (* pi_ str))
  241.       (setf str (/ 1.0 (+ ak dnu)))
  242.       (setf qr (* qr str))
  243.       (setf qi (* qi str))
  244.       (setf str (- (* ckr czr) (* cki czi)))
  245.       (setf rak (/ 1.0 ak))
  246.       (setf cki (* (+ (* ckr czi) (* cki czr)) rak))
  247.       (setf ckr (* str rak))
  248.       (setf s1r (+ (- (* ckr fr) (* cki fi)) s1r))
  249.       (setf s1i (+ (* ckr fi) (* cki fr) s1i))
  250.       (setf str (- pr (* fr ak)))
  251.       (setf sti (- pi_ (* fi ak)))
  252.       (setf s2r (+ (- (* ckr str) (* cki sti)) s2r))
  253.       (setf s2i (+ (* ckr sti) (* cki str) s2i))
  254.       (setf a1 (* a1 t1 rak))
  255.       (setf bk (+ bk ak ak 1.0))
  256.       (setf ak (+ ak 1.0))
  257.       (if (> a1 tol) (go label90))
  258.      label100
  259.       (setf kflag 2)
  260.       (setf a1 (+ fnu 1.0))
  261.       (setf ak (* a1 (abs smur)))
  262.       (if (> ak alim) (setf kflag 3))
  263.       (setf str (f2cl-lib:fref cssr (kflag) ((1 3))))
  264.       (setf p2r (* s2r str))
  265.       (setf p2i (* s2i str))
  266.       (multiple-value-bind
  267.           (var-0 var-1 var-2 var-3 var-4 var-5)
  268.           (zmlt p2r p2i rzr rzi s2r s2i)
  269.         (declare (ignore var-0 var-1 var-2 var-3))
  270.         (setf s2r var-4)
  271.         (setf s2i var-5))
  272.       (setf s1r (* s1r str))
  273.       (setf s1i (* s1i str))
  274.       (if (= koded 1) (go label210))
  275.       (multiple-value-bind
  276.           (var-0 var-1 var-2 var-3)
  277.           (zexp zr zi fr fi)
  278.         (declare (ignore var-0 var-1))
  279.         (setf fr var-2)
  280.         (setf fi var-3))
  281.       (multiple-value-bind
  282.           (var-0 var-1 var-2 var-3 var-4 var-5)
  283.           (zmlt s1r s1i fr fi s1r s1i)
  284.         (declare (ignore var-0 var-1 var-2 var-3))
  285.         (setf s1r var-4)
  286.         (setf s1i var-5))
  287.       (multiple-value-bind
  288.           (var-0 var-1 var-2 var-3 var-4 var-5)
  289.           (zmlt s2r s2i fr fi s2r s2i)
  290.         (declare (ignore var-0 var-1 var-2 var-3))
  291.         (setf s2r var-4)
  292.         (setf s2i var-5))
  293.       (go label210)
  294.      label110
  295.       (multiple-value-bind
  296.           (var-0 var-1 var-2 var-3)
  297.           (zsqrt zr zi str sti)
  298.         (declare (ignore var-0 var-1))
  299.         (setf str var-2)
  300.         (setf sti var-3))
  301.       (multiple-value-bind
  302.           (var-0 var-1 var-2 var-3 var-4 var-5)
  303.           (zdiv rthpi czeroi str sti coefr coefi)
  304.         (declare (ignore var-0 var-1 var-2 var-3))
  305.         (setf coefr var-4)
  306.         (setf coefi var-5))
  307.       (setf kflag 2)
  308.       (if (= koded 2) (go label120))
  309.       (if (> zr alim) (go label290))
  310.       (setf str (* (exp (- zr)) (f2cl-lib:fref cssr (kflag) ((1 3)))))
  311.       (setf sti (* (- str) (sin zi)))
  312.       (setf str (* str (cos zi)))
  313.       (multiple-value-bind
  314.           (var-0 var-1 var-2 var-3 var-4 var-5)
  315.           (zmlt coefr coefi str sti coefr coefi)
  316.         (declare (ignore var-0 var-1 var-2 var-3))
  317.         (setf coefr var-4)
  318.         (setf coefi var-5))
  319.      label120
  320.       (if (= (abs dnu) 0.5) (go label300))
  321.       (setf ak (cos (* dpi dnu)))
  322.       (setf ak (coerce (abs ak) 'double-float))
  323.       (if (= ak czeror) (go label300))
  324.       (setf fhs (coerce (abs (- 0.25 dnu2)) 'double-float))
  325.       (if (= fhs czeror) (go label300))
  326.       (setf t1
  327.               (coerce
  328.                (the f2cl-lib:integer4
  329.                     (f2cl-lib:int-sub (f2cl-lib:i1mach 14) 1))
  330.                'double-float))
  331.       (setf t1 (* t1 (f2cl-lib:d1mach 5) 3.321928094))
  332.       (setf t1 (max t1 12.0))
  333.       (setf t1 (min t1 60.0))
  334.       (setf t2 (- (* tth t1) 6.0))
  335.       (if (/= zr 0.0) (go label130))
  336.       (setf t1 hpi)
  337.       (go label140)
  338.      label130
  339.       (setf t1 (f2cl-lib:datan (/ zi zr)))
  340.       (setf t1 (coerce (abs t1) 'double-float))
  341.      label140
  342.       (if (> t2 caz) (go label170))
  343.       (setf etest (/ ak (* dpi caz tol)))
  344.       (setf fk coner)
  345.       (if (< etest coner) (go label180))
  346.       (setf fks ctwor)
  347.       (setf ckr (+ caz caz ctwor))
  348.       (setf p1r czeror)
  349.       (setf p2r coner)
  350.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  351.                     ((> i kmax) nil)
  352.         (tagbody
  353.           (setf ak (/ fhs fks))
  354.           (setf cbr (/ ckr (+ fk coner)))
  355.           (setf ptr p2r)
  356.           (setf p2r (- (* cbr p2r) (* p1r ak)))
  357.           (setf p1r ptr)
  358.           (setf ckr (+ ckr ctwor))
  359.           (setf fks (+ fks fk fk ctwor))
  360.           (setf fhs (+ fhs fk fk))
  361.           (setf fk (+ fk coner))
  362.           (setf str (* (abs p2r) fk))
  363.           (if (< etest str) (go label160))
  364.          label150))
  365.       (go label310)
  366.      label160
  367.       (setf fk (+ fk (* spi t1 (f2cl-lib:fsqrt (/ t2 caz)))))
  368.       (setf fhs (coerce (abs (- 0.25 dnu2)) 'double-float))
  369.       (go label180)
  370.      label170
  371.       (setf a2 (f2cl-lib:fsqrt caz))
  372.       (setf ak (/ (* fpi ak) (* tol (f2cl-lib:fsqrt a2))))
  373.       (setf aa (/ (* 3.0 t1) (+ 1.0 caz)))
  374.       (setf bb (/ (* 14.7 t1) (+ 28.0 caz)))
  375.       (setf ak
  376.               (/
  377.                (+ (f2cl-lib:flog ak)
  378.                   (/ (* caz (cos aa)) (+ 1.0 (* 0.008 caz))))
  379.                (cos bb)))
  380.       (setf fk (+ (/ (* 0.12125 ak ak) caz) 1.5))
  381.      label180
  382.       (setf k (f2cl-lib:int fk))
  383.       (setf fk (coerce (the f2cl-lib:integer4 k) 'double-float))
  384.       (setf fks (* fk fk))
  385.       (setf p1r czeror)
  386.       (setf p1i czeroi)
  387.       (setf p2r tol)
  388.       (setf p2i czeroi)
  389.       (setf csr p2r)
  390.       (setf csi p2i)
  391.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  392.                     ((> i k) nil)
  393.         (tagbody
  394.           (setf a1 (- fks fk))
  395.           (setf ak (/ (+ fks fk) (+ a1 fhs)))
  396.           (setf rak (/ 2.0 (+ fk coner)))
  397.           (setf cbr (* (+ fk zr) rak))
  398.           (setf cbi (* zi rak))
  399.           (setf ptr p2r)
  400.           (setf pti p2i)
  401.           (setf p2r (* (- (* ptr cbr) (* pti cbi) p1r) ak))
  402.           (setf p2i (* (- (+ (* pti cbr) (* ptr cbi)) p1i) ak))
  403.           (setf p1r ptr)
  404.           (setf p1i pti)
  405.           (setf csr (+ csr p2r))
  406.           (setf csi (+ csi p2i))
  407.           (setf fks (+ (- a1 fk) coner))
  408.           (setf fk (- fk coner))
  409.          label190))
  410.       (setf tm (zabs csr csi))
  411.       (setf ptr (/ 1.0 tm))
  412.       (setf s1r (* p2r ptr))
  413.       (setf s1i (* p2i ptr))
  414.       (setf csr (* csr ptr))
  415.       (setf csi (* (- csi) ptr))
  416.       (multiple-value-bind
  417.           (var-0 var-1 var-2 var-3 var-4 var-5)
  418.           (zmlt coefr coefi s1r s1i str sti)
  419.         (declare (ignore var-0 var-1 var-2 var-3))
  420.         (setf str var-4)
  421.         (setf sti var-5))
  422.       (multiple-value-bind
  423.           (var-0 var-1 var-2 var-3 var-4 var-5)
  424.           (zmlt str sti csr csi s1r s1i)
  425.         (declare (ignore var-0 var-1 var-2 var-3))
  426.         (setf s1r var-4)
  427.         (setf s1i var-5))
  428.       (if (or (> inu 0) (> n 1)) (go label200))
  429.       (setf zdr zr)
  430.       (setf zdi zi)
  431.       (if (= iflag 1) (go label270))
  432.       (go label240)
  433.      label200
  434.       (setf tm (zabs p2r p2i))
  435.       (setf ptr (/ 1.0 tm))
  436.       (setf p1r (* p1r ptr))
  437.       (setf p1i (* p1i ptr))
  438.       (setf p2r (* p2r ptr))
  439.       (setf p2i (* (- p2i) ptr))
  440.       (multiple-value-bind
  441.           (var-0 var-1 var-2 var-3 var-4 var-5)
  442.           (zmlt p1r p1i p2r p2i ptr pti)
  443.         (declare (ignore var-0 var-1 var-2 var-3))
  444.         (setf ptr var-4)
  445.         (setf pti var-5))
  446.       (setf str (- (+ dnu 0.5) ptr))
  447.       (setf sti (- pti))
  448.       (multiple-value-bind
  449.           (var-0 var-1 var-2 var-3 var-4 var-5)
  450.           (zdiv str sti zr zi str sti)
  451.         (declare (ignore var-0 var-1 var-2 var-3))
  452.         (setf str var-4)
  453.         (setf sti var-5))
  454.       (setf str (+ str 1.0))
  455.       (multiple-value-bind
  456.           (var-0 var-1 var-2 var-3 var-4 var-5)
  457.           (zmlt str sti s1r s1i s2r s2i)
  458.         (declare (ignore var-0 var-1 var-2 var-3))
  459.         (setf s2r var-4)
  460.         (setf s2i var-5))
  461.      label210
  462.       (setf str (+ dnu 1.0))
  463.       (setf ckr (* str rzr))
  464.       (setf cki (* str rzi))
  465.       (if (= n 1) (setf inu (f2cl-lib:int-sub inu 1)))
  466.       (if (> inu 0) (go label220))
  467.       (if (> n 1) (go label215))
  468.       (setf s1r s2r)
  469.       (setf s1i s2i)
  470.      label215
  471.       (setf zdr zr)
  472.       (setf zdi zi)
  473.       (if (= iflag 1) (go label270))
  474.       (go label240)
  475.      label220
  476.       (setf inub 1)
  477.       (if (= iflag 1) (go label261))
  478.      label225
  479.       (setf p1r (f2cl-lib:fref csrr (kflag) ((1 3))))
  480.       (setf ascle (f2cl-lib:fref bry (kflag) ((1 3))))
  481.       (f2cl-lib:fdo (i inub (f2cl-lib:int-add i 1))
  482.                     ((> i inu) nil)
  483.         (tagbody
  484.           (setf str s2r)
  485.           (setf sti s2i)
  486.           (setf s2r (+ (- (* ckr str) (* cki sti)) s1r))
  487.           (setf s2i (+ (* ckr sti) (* cki str) s1i))
  488.           (setf s1r str)
  489.           (setf s1i sti)
  490.           (setf ckr (+ ckr rzr))
  491.           (setf cki (+ cki rzi))
  492.           (if (>= kflag 3) (go label230))
  493.           (setf p2r (* s2r p1r))
  494.           (setf p2i (* s2i p1r))
  495.           (setf str (coerce (abs p2r) 'double-float))
  496.           (setf sti (coerce (abs p2i) 'double-float))
  497.           (setf p2m (max str sti))
  498.           (if (<= p2m ascle) (go label230))
  499.           (setf kflag (f2cl-lib:int-add kflag 1))
  500.           (setf ascle (f2cl-lib:fref bry (kflag) ((1 3))))
  501.           (setf s1r (* s1r p1r))
  502.           (setf s1i (* s1i p1r))
  503.           (setf s2r p2r)
  504.           (setf s2i p2i)
  505.           (setf str (f2cl-lib:fref cssr (kflag) ((1 3))))
  506.           (setf s1r (* s1r str))
  507.           (setf s1i (* s1i str))
  508.           (setf s2r (* s2r str))
  509.           (setf s2i (* s2i str))
  510.           (setf p1r (f2cl-lib:fref csrr (kflag) ((1 3))))
  511.          label230))
  512.       (if (/= n 1) (go label240))
  513.       (setf s1r s2r)
  514.       (setf s1i s2i)
  515.      label240
  516.       (setf str (f2cl-lib:fref csrr (kflag) ((1 3))))
  517.       (f2cl-lib:fset (f2cl-lib:fref yr (1) ((1 n))) (* s1r str))
  518.       (f2cl-lib:fset (f2cl-lib:fref yi (1) ((1 n))) (* s1i str))
  519.       (if (= n 1) (go end_label))
  520.       (f2cl-lib:fset (f2cl-lib:fref yr (2) ((1 n))) (* s2r str))
  521.       (f2cl-lib:fset (f2cl-lib:fref yi (2) ((1 n))) (* s2i str))
  522.       (if (= n 2) (go end_label))
  523.       (setf kk 2)
  524.      label250
  525.       (setf kk (f2cl-lib:int-add kk 1))
  526.       (if (> kk n) (go end_label))
  527.       (setf p1r (f2cl-lib:fref csrr (kflag) ((1 3))))
  528.       (setf ascle (f2cl-lib:fref bry (kflag) ((1 3))))
  529.       (f2cl-lib:fdo (i kk (f2cl-lib:int-add i 1))
  530.                     ((> i n) nil)
  531.         (tagbody
  532.           (setf p2r s2r)
  533.           (setf p2i s2i)
  534.           (setf s2r (+ (- (* ckr p2r) (* cki p2i)) s1r))
  535.           (setf s2i (+ (* cki p2r) (* ckr p2i) s1i))
  536.           (setf s1r p2r)
  537.           (setf s1i p2i)
  538.           (setf ckr (+ ckr rzr))
  539.           (setf cki (+ cki rzi))
  540.           (setf p2r (* s2r p1r))
  541.           (setf p2i (* s2i p1r))
  542.           (f2cl-lib:fset (f2cl-lib:fref yr (i) ((1 n))) p2r)
  543.           (f2cl-lib:fset (f2cl-lib:fref yi (i) ((1 n))) p2i)
  544.           (if (>= kflag 3) (go label260))
  545.           (setf str (coerce (abs p2r) 'double-float))
  546.           (setf sti (coerce (abs p2i) 'double-float))
  547.           (setf p2m (max str sti))
  548.           (if (<= p2m ascle) (go label260))
  549.           (setf kflag (f2cl-lib:int-add kflag 1))
  550.           (setf ascle (f2cl-lib:fref bry (kflag) ((1 3))))
  551.           (setf s1r (* s1r p1r))
  552.           (setf s1i (* s1i p1r))
  553.           (setf s2r p2r)
  554.           (setf s2i p2i)
  555.           (setf str (f2cl-lib:fref cssr (kflag) ((1 3))))
  556.           (setf s1r (* s1r str))
  557.           (setf s1i (* s1i str))
  558.           (setf s2r (* s2r str))
  559.           (setf s2i (* s2i str))
  560.           (setf p1r (f2cl-lib:fref csrr (kflag) ((1 3))))
  561.          label260))
  562.       (go end_label)
  563.      label261
  564.       (setf helim (* 0.5 elim))
  565.       (setf elm (exp (- elim)))
  566.       (setf celmr elm)
  567.       (setf ascle (f2cl-lib:fref bry (1) ((1 3))))
  568.       (setf zdr zr)
  569.       (setf zdi zi)
  570.       (setf ic -1)
  571.       (setf j 2)
  572.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  573.                     ((> i inu) nil)
  574.         (tagbody
  575.           (setf str s2r)
  576.           (setf sti s2i)
  577.           (setf s2r (+ (- (* str ckr) (* sti cki)) s1r))
  578.           (setf s2i (+ (* sti ckr) (* str cki) s1i))
  579.           (setf s1r str)
  580.           (setf s1i sti)
  581.           (setf ckr (+ ckr rzr))
  582.           (setf cki (+ cki rzi))
  583.           (setf as (zabs s2r s2i))
  584.           (setf alas (f2cl-lib:flog as))
  585.           (setf p2r (- alas zdr))
  586.           (if (< p2r (- elim)) (go label263))
  587.           (multiple-value-bind
  588.               (var-0 var-1 var-2 var-3 var-4)
  589.               (zlog s2r s2i str sti idum)
  590.             (declare (ignore var-0 var-1))
  591.             (setf str var-2)
  592.             (setf sti var-3)
  593.             (setf idum var-4))
  594.           (setf p2r (- str zdr))
  595.           (setf p2i (- sti zdi))
  596.           (setf p2m (/ (exp p2r) tol))
  597.           (setf p1r (* p2m (cos p2i)))
  598.           (setf p1i (* p2m (sin p2i)))
  599.           (multiple-value-bind
  600.               (var-0 var-1 var-2 var-3 var-4)
  601.               (zuchk p1r p1i nw ascle tol)
  602.             (declare (ignore var-0 var-1 var-3 var-4))
  603.             (setf nw var-2))
  604.           (if (/= nw 0) (go label263))
  605.           (setf j (f2cl-lib:int-sub 3 j))
  606.           (f2cl-lib:fset (f2cl-lib:fref cyr (j) ((1 2))) p1r)
  607.           (f2cl-lib:fset (f2cl-lib:fref cyi (j) ((1 2))) p1i)
  608.           (if (= ic (f2cl-lib:int-sub i 1)) (go label264))
  609.           (setf ic i)
  610.           (go label262)
  611.          label263
  612.           (if (< alas helim) (go label262))
  613.           (setf zdr (- zdr elim))
  614.           (setf s1r (* s1r celmr))
  615.           (setf s1i (* s1i celmr))
  616.           (setf s2r (* s2r celmr))
  617.           (setf s2i (* s2i celmr))
  618.          label262))
  619.       (if (/= n 1) (go label270))
  620.       (setf s1r s2r)
  621.       (setf s1i s2i)
  622.       (go label270)
  623.      label264
  624.       (setf kflag 1)
  625.       (setf inub (f2cl-lib:int-add i 1))
  626.       (setf s2r (f2cl-lib:fref cyr (j) ((1 2))))
  627.       (setf s2i (f2cl-lib:fref cyi (j) ((1 2))))
  628.       (setf j (f2cl-lib:int-sub 3 j))
  629.       (setf s1r (f2cl-lib:fref cyr (j) ((1 2))))
  630.       (setf s1i (f2cl-lib:fref cyi (j) ((1 2))))
  631.       (if (<= inub inu) (go label225))
  632.       (if (/= n 1) (go label240))
  633.       (setf s1r s2r)
  634.       (setf s1i s2i)
  635.       (go label240)
  636.      label270
  637.       (f2cl-lib:fset (f2cl-lib:fref yr (1) ((1 n))) s1r)
  638.       (f2cl-lib:fset (f2cl-lib:fref yi (1) ((1 n))) s1i)
  639.       (if (= n 1) (go label280))
  640.       (f2cl-lib:fset (f2cl-lib:fref yr (2) ((1 n))) s2r)
  641.       (f2cl-lib:fset (f2cl-lib:fref yi (2) ((1 n))) s2i)
  642.      label280
  643.       (setf ascle (f2cl-lib:fref bry (1) ((1 3))))
  644.       (multiple-value-bind
  645.           (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10
  646.            var-11)
  647.           (zkscl zdr zdi fnu n yr yi nz rzr rzi ascle tol elim)
  648.         (declare
  649.          (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-7 var-8 var-9 var-10
  650.           var-11))
  651.         (setf nz var-6))
  652.       (setf inu (f2cl-lib:int-sub n nz))
  653.       (if (<= inu 0) (go end_label))
  654.       (setf kk (f2cl-lib:int-add nz 1))
  655.       (setf s1r (f2cl-lib:fref yr (kk) ((1 n))))
  656.       (setf s1i (f2cl-lib:fref yi (kk) ((1 n))))
  657.       (f2cl-lib:fset (f2cl-lib:fref yr (kk) ((1 n)))
  658.                      (* s1r (f2cl-lib:fref csrr (1) ((1 3)))))
  659.       (f2cl-lib:fset (f2cl-lib:fref yi (kk) ((1 n)))
  660.                      (* s1i (f2cl-lib:fref csrr (1) ((1 3)))))
  661.       (if (= inu 1) (go end_label))
  662.       (setf kk (f2cl-lib:int-add nz 2))
  663.       (setf s2r (f2cl-lib:fref yr (kk) ((1 n))))
  664.       (setf s2i (f2cl-lib:fref yi (kk) ((1 n))))
  665.       (f2cl-lib:fset (f2cl-lib:fref yr (kk) ((1 n)))
  666.                      (* s2r (f2cl-lib:fref csrr (1) ((1 3)))))
  667.       (f2cl-lib:fset (f2cl-lib:fref yi (kk) ((1 n)))
  668.                      (* s2i (f2cl-lib:fref csrr (1) ((1 3)))))
  669.       (if (= inu 2) (go end_label))
  670.       (setf t2 (+ fnu (f2cl-lib:int-sub kk 1)))
  671.       (setf ckr (* t2 rzr))
  672.       (setf cki (* t2 rzi))
  673.       (setf kflag 1)
  674.       (go label250)
  675.      label290
  676.       (setf koded 2)
  677.       (setf iflag 1)
  678.       (setf kflag 2)
  679.       (go label120)
  680.      label300
  681.       (setf s1r coefr)
  682.       (setf s1i coefi)
  683.       (setf s2r coefr)
  684.       (setf s2i coefi)
  685.       (go label210)
  686.      label310
  687.       (setf nz -2)
  688.       (go end_label)
  689.      end_label
  690.       (return (values nil nil nil nil nil nil nil nz nil nil nil)))))
  691.  
  692.